Nonequilibrium density matrix description of steady state 

quantum transport 

Abhishek Dhar * Keiji Saito ^and Peter Hanggi ^ 
January 5, 2012 

Abstract 

With this work we investigate the stationary nonequilibrium density matrix of current 
carrying nonequilibrium steady states of in-between quantum systems that are connected 
to reservoirs. We describe the analytical procedure to obtain the explicit result for the 
reduced density matrix of quantum transport when the system, the connecting reservoirs 
and, as well, the system-reservoir interactions are described by quadratic Hamiltonians. 
Our procedure is detailed for both, electronic transport described by the tight-binding 
Hamiltonian and for phonon transport described by harmonic Hamiltonians. For the 
special case of weak system-reservoir couplings, a more detailed description of the steady- 
state density matrix is obtained. Several paradigm transport setups for inter-electrode 
electron transport and low-dimensional phonon heat flux are elucidated. 
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1 Introduction 



The theory of equilibrium statistical mechanics, as pioneered by Boltzmann and Gibbs, pro- 
vides the prescription for the appropriate density matrix (or density operator)-description of a 
system that is kept under various external constraints. Thus, for systems kept in isolation the 
microcanonical distribution yields the appropriate density matrix, while for systems in weak 
contact with a thermal and particle reservoir the grand-canonical density matrix describes the 
statistical state of the system. For classical systems, equilibrium statistical physics is governed 
by the phase space distribution of the system. A knowledge of the density matrix or the phase 
space distribution then enables one to find various equilibrium and also close to equilibrium 
properties of a system, as exemplified, for example, via linear response theory. 

For systems taken far away from equilibrium there exists no general procedure in obtaining 
its density matrix. Particularly, this holds true for systems that have reached steady states. 
For classical Hamiltonian systems described by a Markovian stochastic dynamics the steady 
state is determined by the stationary solution of the corresponding master equation; e.g. it 
is given by the stationary probability density of the Fokker-Planck generator for continuous 
Markovian processes [1]. Apart from specific situations, however, for example (i) in the presence 
of symmetries such as (strict) detailed balance, or (ii) a single variable state space dynamics 
[1], the task of finding the closed form solution of such master equations presents a profound 
challenge which typically can be obtained only by the usage of extensive numerical simulations 
or algorithms. 

In this context we remind the readers that even for the case of a system being in contact 
with a single bath the corresponding canonical equilibrium is typically not of the common 
Boltzmann-Gibbs structure, as encoded with the exponential of the (negative) bare system 
Hamiltonian and inverse temperature. The latter structure holds rigorously true for weak 
coupling. In presence of strong coupling, however, the corresponding thermal (generalized 
canonical) density operator then typically involves a temperature-dependent "Hamiltonian of 
mean force" [2] which includes entropic contributions that explicitly depend on the system-bath 
coupling strength. 

Regrettably, no such general concept as the canonical Boltzmann-Gibbs density matrix 
structure in terms of the bare (or even modified Hamiltonian of mean force) is available when 
the open system is subjected to steady state transport. Put differently, there are no generic 
results known for stationary nonequilibrium statistics. This latter situation in fact is not only 
substantially more complex but presently is also less researched. It is thus of outmost impor- 
tance to gain further insight into this objective of obtaining the underlying nonequilibrium 
density matrices that govern quantum and/or classical transport. For example, the explicit 
form of a corresponding nonequilibrium density matrix not only determines the linear response 
due to an additional external perturbation of such a nonequilibrium steady state (NESS), but 
also its higher order response functions. 

A particular, exactly solvable case is that of heat conduction occurring in a one- dimensional 
ordered harmonic chain when connected to two baths at different temperatures. If the two 
baths are modeled therein as being stochastic with corresponding stochastic forces acting on 
the system of interest, the exact nonequilibrium steady state phase space distribution for 
this problem was evaluated by Rieder, Lebowitz and Lieb [3]. An extension to the case of 
higher dimensions was later obtained by Nakazawa [4]. Heat conduction in quantum harmonic 
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oscillator chains has been studied by several authors [5, 6, 7, 8], but thus far no explicit results 
are known for the precise form of the quantum mechanical steady state density matrix. Some 
formal results for the NESS density matrix of general quantum mechanical systems have been 
obtained in the works of Zubarev [9] and McLennan [10] and have more recently been discussed 
in specific models [11, 12, 13]. 

The problem of obtaining the NESS in explicit form presents a formidable challenge already 
for classical open systems (see above discussion). This objective therefore is typically even more 
intricate for a quantum NESS. Indeed, in presence of general nonlinear interactions this task 
is simply inaccessible without invoking also extensive and cumbersome numerical means and 
methods. To obtain general analytic insight over whole parameter ranges thus necessitate to 
confine the objective to stylized situations only, that allow for explicit closed form calculations. 

With this work, we consider generic setups for steady state quantum transport described by 
a bilinear Hamiltonian. The aim is to find systematic procedure for obtaining explicit results 
for the NESS density matrix for this class of systems. We demonstrate that it is possible to 
obtain the complete NESS density matrix explicitly. We also show that when the coupling 
strength between the system and reservoirs are extremely weak, the NESS density matrix is 
given by an effective Gibbs state where each mode is formally only in equilibrum with a mode- 
dependent effective temperature which depends, however, in a complex manner on both bath 
temperatures. 

We consider two generic setups for stationary nonequilibrium quantum transport, a first 
one involving fermions and the other one bosons as carriers. The first setup consists of elec- 
tron and heat transport in a fermion setup of non-interacting particles which are connected 
to fermionic baths at different temperatures and chemical potentials. The second scenario 
consists of heat conduction occurring in harmonic crystals connected to oscillator baths kept 
at different temperatures. For both these problems it is known from prior studies, using vari- 
ous approaches such as the nonequilibrium Green's function formalism [14, 15], the quantum 
Langevin equations approach [16, 17, 18], and the C*-algebra approach [19], that it is possible 
to express all two point correlations in the NESS in terms of appropriate Green's functions. 
Because these systems are non-interacting it is evident that the two-point correlations contain 
necessary information on all higher-point correlations and hence should completely specify the 
NESS. With this study we give the procedure for finding the explicit NESS density matrices 
from a knowledge of the two-point correlations in these systems. 

For the case of a weak-coupling among system and the baths we are able to obtain explicit 
results. We further present explicit examples in simple one-dimensional models which illustrate 
our general procedure and also demonstrate the accuracy of the weak-coupling approximation. 

The outline of the paper is as follows. In Sec. (2) we present the general procedure for 
construction of the NESS for the electron and phonon transport problems. The special case of 
weak coupling between system-bath is discussed in Sec. (3). In Sec. (4) we discuss some illus- 
trative examples of models where both system and reservoirs are taken to be one-dimensional 
chains. Finally we end with a discussion in Sec. (5). 

2 Construction of Steady State Density Matrix 

In this section, we outline the general procedure to obtain the steady state density matrix in 
quantum transport described by a bilinear Hamiltonian. We focus on electric conduction as an 
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example of fermionic transport, and phononic heat conduction as bosonic transport. Because 
the resulting NESS density matrix becomes be Gaussian, the pertinent procedure fist is in 
finding the explicit form of the two-point correlation functions of physical quantities and next 
relating these to the Gaussian distribution. 



2.1 Steady state density matrix for non-interacting electron trans- 
port 

We consider the typical setup of transport in the Landauer approach wherein a system is 
connected to two reservoirs initially kept at different temperatures and chemical potentials. 
At long times the system reaches a nonequilibrium steady state with a mean rate of flow of 
charge and energy current. One starts out by writing the full Hamiltonian of the system plus 
reservoirs and here we consider a tight-binding approach of non-interacting electrons. We use 
the following notation: for sites on the system (S) we shall use the integer indices l,m,n, ■ ■ • ; 
for sites on the left reservoir (L) we employ the Greek indices a, v\ and finally, for sites on 
the right reservoir (R) we use the primed Greek indices a', v' . We consider quantum transport 
with the following overall Hamiltonian reading: 

T-L = Us + Hl + Hr + Hls + Urs , (2.1) 
where H s = ^ H i™ 4 C ™, Ml = ^2 H ^ C " c <" Ur = S C " ,cv ' 

lm av ot'v 1 



la 

H RS = J2 H ™ t ^a' + Old } , 



la' 

where c\ c denote creation and annihilation operators satisfying the usual fermionic anti- 
commutation rules and we assume that the matrices H, H L , H R are symmetric and real-valued 
while H LS , H RS are real-valued. In the above setup we assume that the system possesses a 
finite number of lattice sites iV while the left and right reservoirs have and Nr sites which 
will eventually be made infinite. The parts Hs, Hl and Hr denote the Hamiltonians of the 
isolated system, left and right reservoirs respectively, while %l and 1-Lr describe the coupling 
of the left and right reservoirs to the system, which have been taken to be real. To obtain a 
NESS for the system we consider an initial state at time t = t given by the following product 
density matrix: 

p(t )=Ps®P°L®pR, (2-2) 

where p\ ~ e~^ i ~ Mi - v ' i ^ Ti (p° R ~ e ~(' H R-i 1 R J ^R)/ T R^ [ s the equilibrium grand-canonical density 
matrix of the left (right) reservoir, corresponding to temperature Tl (Tr) and chemical po- 
tential pi (pr) with Nl,Nr the total number operators, while p% denotes an arbitrary initial 
density matrix for the system. We then time-evolve the whole system with the full Hamiltonian 
given in Eq. (2.1) so that at time t the full density matrix is given by: 

p{t) = eint-toyhp^-mit-toyn (23) 



4 



Our principal objective is the long time limit of the steady state reduced density matrix for 
the system under consideration, i.e., 

p s = lim lim Tr L>R p(t) , (2.4) 

to— oo baths— »oo 

where the trace, Tr l,r, is over all the degrees of freedom of the two baths. In doing so we 
implicitly assume that the quantum transport setup is so that (i) it possesses a long time 
limit in the limit of infinite many bath degrees of freedom and (ii) that the interactions within 
the quantum system and the interaction with the bath degrees of freedom are such that the 
emerging asymptotic nonequilibrium steady state density matrix indeed is time- independent. 
Let us introduce the two-point correlation function 

(cIc 1 ) = Tt s [cIc iPs ] = Tt[cIc iP } (2.5) 

where the first trace, Tis, is over system degrees of freedom and the second trace is over all 
degrees of freedom. Because of the quadratic form of the total Hamiltonian the two-point 
correlations in the NESS can be exactly calculated using various methods [20, 21, 16]. 
The correlations can be expressed in terms of the following Green function: 

where and are self-energy terms which model the effect of the infinite reservoirs on 
the isolated system Hamiltonian. The self energies can be written in terms of the isolated 
reservoir Green functions g^{oS) = [fku ± ie — H 1 ]' 1 , g%{u) = [fvjj ± ie — H 11 ]^ 1 and the 
coupling matrices H LS and H RS , reading 

E±H = H LS g±(to)H LS \ E±H = H RS g±(u;)H RS ' . (2.7) 

Let us next define T L (u) = [S~ -S+ }/{2i) , T R (u) = [E R - £+]/(2i). With these definitions 
one finds the following expressions for the steady state correlation matrix: 

c ml = ( 4, cj ) 

= / dw - [ (G + T L G') lm f(u,n L ,T L ) + (G + T R G~)i m f(u,fi R ,T R ) ] , (2.8) 

where f(u,(j,,T) = l/fe' 3 ^ - ^ + 1] denotes the Fermi function. 

We demonstrate next how the NESS density matrix ps can be fully expressed in terms of 
these correlations. Note that the matrix C is Hermitian, since at any time Tr[ c\ c m p(t)] = 
Tr[ c^j q p(t) ]*, where (*) indicates complex conjugation. This result can also be directly 
verified from the form in Eq. (2.8). Consequently the matrix C can be diagonalized with a 
unitary matrix U to read: 

U ] CU = D = Diag(di,d 2 ,--- ,d N - U d N ), (2.9) 

where the matrix D is the diagonal matrix. Using the unitary transformation, we define new 
fermionic operators as 

c' s = ]T[/ M q, s=l,...,N. (2.10) 
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Obviously, these new fermionic operators preserve the anti-commutation relations, {c' s ,c'J,} = 
5 StS >. The steady state density matrix is a diagonal matrix in terms of these new fermion 
operators. Note that the two-point correlation of new fermionic operators read {d]d s ,) = 5 SyS 'd s . 
From this we find the corresponding effective Fermi-Dirac distribution for each fermion s. 
Consequently, the steady state matrix p s is formally given by 



Ps 



nexp \— a s c'VJ 



A [1 + eX p(-a s )] 

exp j^— X^/,rrt C lAi >m C n 



, (2.12) 
IL=i t 1 + ex P(- fl s )] 

A = J7* Diag(ai, a 2 , • • • ,a N _i,a N )U T , (2.13) 

a s = In (d; 1 - 1) . (2.14) 

To obtain Eq. (2.14) we used the relation ( c'} c' s ) = d s = l/[exp(a s ) + 1]. This completes 
our derivation of the expression for the steady state density matrix for noninteracting electron 
transport. 



2.2 Steady state density matrix for noninteracting phonon trans- 
port 

We next consider heat conduction in general harmonic networks. Examples of such a system 
are dielectric crystals for which the harmonic crystal provides a very good description. As 
before we again consider the usual Landauer-like framework of a system connected to two 
reservoirs kept at different temperatures [7, 8]. The reservoirs are themselves modeled as 
collections of harmonic oscillators. Let us assume that the system has iV Cartesian positional 
degrees of freedom {xi}, I = 1,2 ... ,N with corresponding momenta {pi}. These satisfy the 
usual commutation relations [xi,p m ] = ih5^ m and [x/,x m ] = [pi,p m ] = 0. Similarly the left 
reservoir degrees of freedom are denoted by {x^,p^}, a — 1, ... , N L and the right reservoirs 
by { x a>->Pa'}i ot? = 1, . . . , N R . We will use the vector notation X T = (xi, x 2 , ■ ■ ■ , x^), P t = 
(pi,P2, ■ ■ ■ ,Pn) and similarly X L , X R , P L , P R . The general Hamiltonian for the system coupled 
to harmonic reservoirs is then given by: 

U = H s + n L + H R + H LS + H RS , (2.15) 
where Hs = ^P T M' 1 P + ^X T K X , 

U L = -[P L ] T \M L }- 1 P L + -\X L ] T K L X L , 

H R = l -[P R ] T [MY 1 P R + \[X R ] T K R X R , 
Uls = X T K LS X L , U RS = X T K RS X R , 

where M, M L , M R and K, K L , K R denote respectively the mass matrix and the force- 
constant matrix of the system, left reservoir and right reservoir, while K LS and K RS denote 
the linear coupling coefficients between the two reservoirs and the system. 
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Again we consider the time evolution of the coupled system plus reservoirs starting from 
an initial product density matrix of the form Eq. (2.2) with p\ ~ exp(— 'H.lI^bPl) and p R ~ 
exp(— "H^/ZcbTr) and the system being in an arbitrary initial state. At long times the system 
reaches a NESS described the reduced density matrix ps = Tr^^p(t — >■ oo). In order to 
construct ps , we start with defining the appropriate correlation matrix as in the previous 
section for electron transport. In doing so we consider the 2N x 2N covariance matrix defined 
with the column vector ip = (x±, • • • , Xn,Pi, • • • ,Pn) T '- 



cptp 1 ) =Tr s [ tptp 1 Ps ] . (2.16) 
For this covariance matrix, we write the symmetric and anti-symmetric parts as 

Cs = \{C + C T ) , (2.17) 

Ca = \(C-C T ) = l p (2.18) 

(.: s). w 

where 1 and are respectively N x N identity and zero matrices. The matrix expression 
of anti-symmetric part Ca is automatically determined by commutation relations between 
coordinate and momentum variables. 

The symmetric part of covariance matrix Cs is given by 

c - ( {xxT) 1(xp t + Ipx t ] t }\ r22m 

° 5 - {±(xp t +[px t ] t ) (pp t ) ) ' 1 j 

As for the electron case these correlations are known in terms of the following phonon Green 
function: 

G± = -M^ + K-Ef-E*' (2 ' 21) 

where the self-energies can be expressed in terms of the isolated reservoir Green functions 
gf(u) = [ -M L {u± it) 2 + K L J" 1 , g±(u>) = [ - M R (co ± ief + K R ] _1 and the coupling 
elements K LS , K RS . These self energies thus read 

E±H = K LS gt(u) [K LS ] T , EJ( W ) = K RS g±(u) [K RS ] T . (2.22) 
Defining T L (u) = Im[ S+ ] , T R (u) = Im[ £+ ], we find [17, 22, 23]: 

/oo t) 
duj^- G + T a G- g(u,T a ) , 

■°° a=L,R 

/oo fc 2 

^ — ]T MG T a G M g(u,T a ) , 

i(XP r + [FIT) = r dw— E G + T a G-M «?( W ,T a ) , (2.23) 

2 J-oo K ±f n 



a=L,R 



where g(u,T) = coth(ftw/2k B T). 

We next show how the steady state density matrix can be expressed in terms of the cor- 
relation matrix. For this it is necessary to consider symplectic transformations [24]. We first 
introduce the symplectic matrix S, satisfying 

SJS T = J, (2.24) 
SC S S T = D = Diag (dx,.-. ,d N ,d ir -- ,d N ). (2.25) 

The procedure to find S is detailed in the Appendix (A). 

By using the symplectic transformation with the matrix S, the new operators ip' = 
(x[, ■ ■ ■ , x' N ,p[, ■ ■ ■ ,p' N ) T are defined as: 

N 

<p' s = ^Ss^tpt ,s = l,...,N. (2.26) 
i=i 

The most important property of the the symplectic transformation, following from Eq.(2.24), 
is that it preserves the commutation relations and we have [x S) iv] = ih5 SyS > and [x s ,av] = 
[p s ,Ps>] = 0. The steady state density matrix can then be written in terms of these new 
operators and we end up with the general main result: 



Ps = 



"Q exp[-a s {x' s 2 +p' s 2 )} ^ ^ 

exp [— Lp T Alp\ 



nil z s 



(2.28) 



A = S T Diag (ai,--- ,a N ,a ir -- ,a N ) 5, (2.29) 
Z s = [ 2smh(ha s ) ] _1 , (2.30) 
a s = fT 1 cothr 1 {2d s /h) . (2.31) 

In computing the normalization factor Z s , we have used the second quantization representation 

x' s = \J\{b\ + b s ) ,p' s = — b s ), where b s and b\ satisfy [& s ,feg/] = 5 SjS >. Then we obtain 

the expression a s (x' s 2 + p' s 2 ) = 2ha s (b\b s + 1/2). The relation between d s and a s is then found 
by looking at the averages ( x' 2 ) and ( p' 2 ): 

( x? ) = ( V? ) = \ coth(ha s ) = d s . (2.32) 

Finally we also consider here the classical limit H — > 0. In this limit, we have the simple 
relation d s = l/(2a s ). Then, the matrix A is given by 

A = ^SFD^S (h->0). (2.33) 

Hence, from the relation D 1 = (SCS 7 ^)^ 1 = (S T )~ 1 C~ 1 S~ 1 , we find the following expression 
of the matrix A in the classical limit: 

A = \c- x (h^O). (2.34) 
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Thus we recover the form that is expected for a general Gaussian probability measure. We note 
that in the classical case, for Gaussian white noise reservoirs, the correlation matrix C can be 
explicitly determined for ordered harmonic lattices [3, 4]. For arbitrary harmonic networks, 
they are given by the high temperature limit of Eqs.(2.23), with appropriate choices of the 
bath spectral functions. Finding the inverse of the correlation matrix presents, however, a 
more difficult task. 



3 Weak coupling limit 

In this section, we consider the special case of a weak coupling between the system and reser- 
voirs. We note that it is essential that the weak-coupling limit is taken after the coupled 
system-reservoirs have evolved for an infinte time and thus reached the NESS. Generally, when 
the coupling strength is weak, the density matrix can be expanded in terms of the coupling 
strength. In this case, the zeroth order term in the coupling strength determines the overall 
structure of the electron density profile in the electron conduction case, and the temperature 
profile in the case of phonon heat conduction. The higher order terms of the expansion deter- 
mine the amount of current flowing in the system. Therefore, although the coupling strengths 
must be finite for finite current, even the zeroth order contribution in the expansion of the 
density matrix carries important information on the steady state. In this section, we focus on 
the 0-th order contribution in the weak coupling expansion of the steady state density matrix, 
which we here refer to as the density matrix in the weak coupling limit. We emphasize that at 
no instant we switch off the coupling strength which is always kept finite, but small. 

On decreasing the coupling strength, the current decreases; however, even in the limit of 
zero current, the steady state density matrix is non-trivial and different from the equilibrium 
density matrix. In fact, we will find that the NESS is non-unique in the sense that it depends on 
the way the system-coupling strengths are made to vanish. For the case where the temperatures 
and chemical potentials of the two reservoirs are chosen equal one obtains, in the weak-coupling 
limit, the expected equilibrium grand-canonical (for electron case) and canonical (for phonon 
case) distributions. 



3.1 Electron transport 

We first note that the system's Hermitian Hamiltonian matrix H has the eigenvalue-equation 
Y^m Hi t mV m (s) = X s Vi(s), hence can be diagonalized by the unitary transformation V as 

V j HV = X, V ] V = I. (3.1) 

Next we use the spectral decomposition: 

g+ = yy- 1 ^ -jf -£+-£+] - 1 \v ] ]- l v ] 

= V[hu- X- Vt(E+ + E+)V]" 1 V j . (3.2) 

From this it follows that in the weak coupling limit ££, — > 0, the matrix element G^ m is 
effectively given by 

I,m ~ ^ hu + x s - i( s | r | s ) ' 1 j 
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where ( s \ T \ s' ) = *Ei,m V i*( s )( T )i,mV m (s') and T = T L + T R . It can be shown that the 
off-diagonal terms of the inverse matrix in Eq. (3.2) are of the order of the coupling strength. 
This contribution disappears, however, in the following calculation of the correlation function, 
given this weak coupling limit. The real part of is negligible compared to the remaining 
real parts and thus can be dropped. Hence we obtain: 

/OO fc 
00 a=L,R j,k 

J U n^ ^ hu-X s - i (s\r\s) [Ta] ^ 

00 a=L,Rs,s',j,k s 

nw — A s r + i{ s' I 1 I s' ) 

A careful examination of the limit ( s \ F a \ s ) — > exhibits that only the terms s = s' survive 
in the above summation, yielding: 

/„t„\ r 1, n YM ( f I EM I s > rt, ■ ■ t \ 

{ cl c, > = y i-LL - + < « I rM I . >» T " ] ■ 

Oh — Li, IX S 

Next, making use of the identity 



we find: 



lim = n5(x - a) , 

e^o ( x - a) 2 + e 2 



(4q) = ^^( S )K( S )e s 

s 

where e s = ^ y^t^tjj f(^s/h, /2 a ,T a ) 

a=L,R \ I I / 
= 1L f{\ s /ti,HL,T L ) +7fl f(X s /h,/J, R ,T R ) , 

( g I r L 1 s ) ( s 1 r fl 1 s ) 

where 7 L = - — — — — - , 7^ = - — — — — - = 1 - 7 L . 



Note that, in the above expression, the limit ( s \ T a \ s ) — > is always implied and it is 
then evident that the ratios 7l,7r depend on the way the couplings — > 0. From the form 
above we can interpret e s as an effective occupation probability of the energy-level X s of the 
isolated system and this probability depends on the temperatures and chemical potentials of 
the two reservoirs. Defining the diagonal matrix E with elements e s , we have V^CV = E . 
Comparing with Eq. (2.9) we see that the same unitary transformation which diagonalizes H 
also diagonalizes the correlation matrix C and we have U = V, D = E. 

Using the results in Eqs. (2.13, 2.14) we then find, a s = ^(e^ 1 — 1), and A = V* Diag (ai,a 2 , . . . , a N ) V 1 
which in turn yields the steady state density matrix in Eq. (2.13). For the equilibrium case 

= Hr = T L = T R = T we have d s = e s = f(X s ,/i,T) , hence a s = (A s - /i)/(k B T) and 
A = [H — fiI]/(kBT), as expected. Thus we obtain, the non-trivial result, that the density 
matrix of a system, weakly coupled to two reservoirs at the same temperatures and chemical 
potentials, is given by the grand-canonical distribution of the isolated system. Note that this 
is not the case for the case of strong coupling. 
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3.2 Phonon transport 



For the harmonic model we first note that there exists a real normal mode transformation 
matrix V, with elements V/(s)which satisfies: 



V T MV 



v t kv = n z , 



where f2 is the diagonal matrix with elements as normal mode frequencies. It is easily verified 
that the matrix 



S = 



_n-i/2y T 

n l ' 2 v T M o 



is symplectic i.e SS T = J and further has the following property: 

K 1 



M 



n- 1 o 
o n- 1 



(3.4) 



(3.5) 



We now show that the correlations for the harmonic system in the weak coupling limit are 
given by: 



(x x T ) = v n- 1 ' 2 e n- 1 ' 2 v T 
(X p T + [P x T f) = o 

(P P T ) = M V fi 1/2 E fi 1/2 V T M, 
where we have defined the diagonal matrix E whose elements are given by: 



e s = 



h \ ^ 

2 ^ 

a=L,R 
ft 

7 L - coth 





T a 


IO 




r | 


s) 



coth 



2k B T a 



m 8 

2k B T L 



ft 

+ 1r g coth 



where = 







IO 




1 r | 


O 



1R = 



(0 




IO 


(« 


1 r | 


o 



s=l,2,...,iV 
2k B T R _ ' 

1 - lL ■ 



(3.6) 
(3.7) 
(3.8) 



(3.9) 



Let us define the effective temperature T s for each normal mode through the relation, reading: 



h 



e, = 



1 



2 
2k B 



coth 



l2k B f s 



coth 



7l coth 



giving 



(3.10) 



2k B T L 



+ 7i? coth 



2k B T R 



which notably depends on both temperatures Tl and T B . We remark here again that in the 
above expressions the limit ( s \ T a \ s ) — > is implied and it is then clear that the ratios jl, 
7b depend on the way the couplings — > 0. 

To prove the above results, Eqs. (3.6, 3.7, 3.8, 3.9), we first introduce the following spectral 
decomposition: 

G+H = VV- l [-Mu 2 + K -T l + L -T l + R ]- l [V T ]- l V T 

= V [V T ( - Mu 2 + K - £+ - £+ ) V V T 

= V [ -co 2 + n 2 - v t -e+v- v^+v}- 1 V T . 
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From this it follows that in the weak coupling limit ££, Sjj — y 0, the matrix element G+ m is 
effectively given by 

r+ - V Vi{s)v m {s) , v 

S-L_ u 2 + fi 2_^ | r | s ) ' 1 j 

s 

where ( s | T | s' ) = S/m^( s ) I\ m V rn (s') and T = T L + r R . It can be shown that the 
off-diagonal terms ( s | T \ s' ) for s ^ s', as well as the real part of give lower order 

contributions in the weak coupling limit and can be dropped. Hence we find: 



■oo Z7F r d „■ u 



00 " a=L,R j,k 



X ^)K^) , y , 

A careful examination of the limit ( s | r a | s ) — y shows that only the terms s = s' in the 
above summation survive and we then obtain: 

[~h, TT Vjjs) ( s I T a (u) I s ) V m (s) 
{ 1 m) J 00 2n (-^ + tt 2 ) 2 + ( s I r(w) I S > 2 Ia) ■ 

Now we note the following identity: 

lim — 1 = ^- [5(x - a) + 5(x + a)} . (3.12) 

e^o (x 2 - a 2 ) 2 + e 2 2a 

Using this and the fact that T a (u) and g(u>) are both odd functions of u, one arrives at: 

( *« x m ) = ^ h 2 Vi{s)v m {s) £ ( /Jiri' s g ^ ( ^' Ta) ' (3 - 13) 

s a=L,R \ I I / s 

which proves Eq. (3.6). Similarly we can evaluate other correlations and obtain Eqs. (3.7,3.8). 

From the form of the correlations in Eqs. (3.6, 3.7, 3.8) we deduce that the matrix S 
given in Eq. (3.4) provides the required symplectic transformation in Eq. (2.25) with D = E. 
Therefore, using Eq. (2.31) and the definition in Eq. (3.10) we get a s = tt s /(2kBT s ). Finally 
Eq. (2.29) gives A = S T nf- 1 S/(2k B ) and then from Eq. (2.28) we obtain p s . This density 
matrix corresponds to each of the normal modes of the harmonic system being in equilibrium 
at an effective temperature T s . For the equilibrium case Tl = Tr = T, we find, using Eq. (3.5), 
(p T A(p = V, s /(k B T). This result is expected but is non-trivial, and is valid only in the weak- 
coupling limit. 

4 Application to generic setups 

4.1 Electron transport in a one-dimensional wire 
4.1.1 System with single site 

We consider the system plus reservoir to consist of a single site, such as e.g. realized with 
a single-level quantum dot, that is connected to two one-dimensional reservoirs. The full 
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Hamiltonian then reads: 



% — Us + %l + Hr + H-ls + Urs j 
where U s = ecjco , 

oo oo 

Hi = ~J2t[ 4ca+i + 4 +1 Ca ] , -Hi? = - ^ *[ C L' C «'+1 + C L+1 C «' ] > 
a=l a'=l 

"Hl5 = c Ll c + 4 C "=1 ] , ^i?5 = ~4t C 1'=1 C + 4 C a'=l ] ■ I 4 - 1 ) 

The self-energies can be expressed in terms of the Green functions of the uncoupled reservoir 
Hamiltonian g\ R and the coupling elements t' LR . Defining u = — 2tcosg, where < q < ir, 
we find that for \u\ < 2t: 

= -X*, E+(o/) = -^e*, (4.2) 

r+H = ^sing, r+H = ^sing. 
Hence the system's Green function emerges to read: 

G » = ^- e -E>)-E^) ' 

The correlation matrix element for the single-site problem is then readily obtained, reading 
given by: 

f 2t h 

d=(clc }= / dw-\G+(u)\ 2 [V L {u) f(oo,fi L ,T L ) + T R (co) f(u,fi R ,T R ) } . (4.4) 

J-2t 71 

Consequently we find for the steady state nonequilibrium density matrix for this case the 
explicit result 



exp(— ac ^c ) 



Ps = Z U (4.5) 

1 + exp(— a) 

where a = ln(d _1 — 1) . 



4.1.2 System composed of two sites 

We next consider a system where the reservoirs are identical to those in the previous section, 
while the system Hamiltonian and system-bath couplings are as follows: 

Us = ei4 C l + e 2<4 C 2 - t( C l C 2 + 4 C l) 

Uls = -t' L [ 4=i c i + c\ c a=i ] , Urs = -t' R [ cl, =1 c 2 + c{c a > = i ] . (4.6) 
The self energies are again given by Eq. (4.2) and the system's Green function is then 
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In this case it is difficult to construct explicitly the required unitary matrix U though it is 
straight-forward to evaluate it numerically and from that find the steady state density matrix 
given by Eq. (2.11). 

In what follows we present numerical precise results for this setup. In our numerics we 
use the following set of parameter values: t = l.0,t' L = t' R = 0.05, ei = 0.2, e 2 = 0.4, T L = 
0.25, T R = 0.25. The right reservoir chemical potential is fixed at p: R = 0.0 and we study the 
NESS for different values of A/x = /x^ — /x#. 

The Green function in Eq. (4.7) is first obtained and then all the elements of the correlation 
matrix given by Eqs. (2.8) are evaluated by numerical integration. As examples we give below 
the correlation matrices for the equilibrium case A/i = and for A/x = 2.0. 

„ / 0.519 0.465 \ 

Cs = { 0.465 0.427 ) = ' 

„ / 0.726 0.271 + xU000473 \ „ . 

Cs ~ { 0.271 - x0.000473 0.672 ) f ° r ^ ~ 2 '° " 

The electron current in the chain is given by j e = 2t/m[(c{c 2 }] and in the above example 
j e = 0.000946. 

As discussed in Sec. (2.1) the NESS density matrix assumes the form: 

exp(— Ac) 

PS = [l + exp(-a 1 )] [l+exp(-a 2 )] ' ( " 8) 

where c = (ci,c 2 ) T and we numerically determined the coefficients ai,a 2 and the matrix A. 
Finding the eigenvalues and eigenvectors of C yields the matrix D and the unitary matrix 
U, respectively. We evaluate a\ = ln^ 1 — 1), a 2 = ln(<i 2 1 — 1) and numerically obtain the 
steady state matrix 

A = U* Diag (oi,a 2 ) U T . 

Note that for A/x = (/x^ = /ir = 0) and with a weak-coupling to reservoirs, we expect 
the result, ps = p eq ~ e -PCHs-^) an( j h ence 

0.8 -4.0 
-4.0 1.6 



A eq 



In Fig. (1) we depict the matrix elements An, A 22 and i?e[Ai 2 ] as functions of the chemical 
potential difference A/x. In the inset we also evaluated the electron current; i.e. j e = 2 I m[C 12] 
and show as well im[A 12 ] . 

The matrix elements An, A 22 and the real part of A 12 can be obtained from our analyti- 
cal weak-coupling results in Sec. (3.1). First we obtain the eigenvalues X s and eigenfunctions 
Vi(s),s = 1,2, corresponding to the isolated system Hamiltonian Us- This provides the re- 
quired unitary transformation which diagonalises the matrix C. For the present two-site setup 
the corresponding eigenvalues, which determine the matrix elements of D, generally given by 
Eqs. (3.4), take on the following form: 



j l 1 y n j 7i ^ 1 

(J'S . 9, t \ 1 r\ . O 1 , n ,„ (\ .. \ Irr * ' 





vAs)? 1 


f 2 


vM\ 2 + t> R 2 \ 


V 2 (s)\ 2 e^-^)/T L + 1 







f 2 







for s = 1, 2. The weak-coupling results for An, A 22 and i?e[Ai 2 ] are depicted in Fig. (1) with 
dashed lines. We notice that these are is excellent agreement with the values obtained from 
exact numerics. 
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Figure 1: (color online). Plot of the NESS matrix elements A as a function of the chemical 
potential difference A/x = Hl — with fixed [Mr = 0.0 and with the remaining parame- 
ters as given in the text. The dashed lines depict results obtained from the weak-coupling 
approximation. The inset shows the electron current j e = 2Im[Ci2\ together with im [A12]. 
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4.2 Phonon transport in one-dimensional oscillator chain 
4.2.1 System consisting of a single oscillator 

We consider our system plus reservoir to be described by the full Hamiltonian 



2M 2 

N 



Pa i k(x a iCa+l) 2 ! k L {x a= \ x) 2 

1 ' 

N 



2m 

a=l 



p 2 a , k(x a , - x a , + if k' R (x a ,=i - xf 



9ttj 9 



2m 

a'=l 

where we assume x a= N+i = x a i=N+i = 0. The above Hamiltonian can be written in the 
canonical form: 

n = Us + Ul + Ur + Uls + Urs , (4.9) 

P 2 (K + K + k' R )x 2 
where n s = ^ + < ^ • 

^ " 2^ + 2 + ' 

Q=l 

a/ Pa' , k(x a > — x a > + i) 2 k' R x 2 a , =l 
n « - 2^ + 2 + ' 

a'=l 

"Hls = -k' L x a=1 x, U RS = -k' R x a , =1 x . (4.10) 

The self-energies can be expressed in terms of the Green functions of the uncoupled reservoir 
Hamiltonian g^ R (oj) and the coupling elements k' LR . We define u 2 = (2k/ m) (1 — cos 5), where 
< q < 7r. Then, we find that for \u\ < u m = 2(k/m) 1 / 2 : 

, k' L 2 cos g - (1 - u L ) + i sin g , _ A:^ 2 cos q - (1 - + i sin g 

k 2(1 -u L )(l- cos q)+u\ ' ^ = fc 2(l-u fl )(l-cosg)+«| l ' 

r ^ ^! ^ r y = tf ^9 

MJ fc 2(l- ML )(l-cosg)+ M | ' M; k 2(1 -u R )(l- cos q)+u 2 R ' 

(4.11) 

where ml = fc^/fc and = Hence the Green function is given by: 

1 

-Moo 2 + k Q + k' L + k' R - E£(w) - £+H 



<?» = ^ ; ; TJ ■ ^ ^ • ( 4 - 12 ) 



It is not difficult to verify that T(ui) = 4T l(uj)T R (u)\G + (oj)\ 2 gives the correct transmission 
coefficient as can be independently obtained by evaluating the transmission of plane waves 
from the left reservoir to the right one, across the intermediate system. 
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The correlation matrix elements for the single-particle problem are obtained as: 
ci = (x 2 ) = I m du- \G + {u)\ 2 [T L (u) g(u,T L ) + T R (u) g(u,T R )] , 



o 



7T 



C2 = (P 2 ) = / du \G + (lj)\ 2 [T L (u) g(u,T L ) + T r (uj) g(uj,T R ) ] , 

Jo 71 
(xp + px) = , 

where co m = 2{k/m) 1 / 2 and g(tu,T) = coth([3hu/2). Using the prescription in Sec. (2.2) we 
find that d\ = (cic^) 1 / 2 and 

S= ( ,. ,J - <Cl/C2)1/ I ) • (4-13) 



yielding the explicit NESS density matrix: 

-[A 11 x 2 +A 22 p 2 ] 



z 

1/2 / \ 1/2 

C 2 \ , I C l 



PS = 

where An = [ — I a , A 22 = I — ] a 
\ c iJ \ C 2J 

a = h- 1 coth- 1 ^" 1 ^!^) 1 / 2 ] , 
Z = psinh^a)]- 1 . 

4.2.2 System composed of two coupled oscillators 

In this case the baths have the same Hamiltonians as in the previous section while the system 
Hamiltonian and system-bath couplings are given by: 

_Pi_ j|_ (h + k' L )x 2 k{x l -x 2 f (fc 2 + fc^x 2 , 
5 2mi 2m 2 2 2 2 

= -k' L x a=1 x 1 , U RS = -k R x a > =1 x 2 ■ (4.14) 

The self-energies are again given by Eq. (4.11) and the system's Green function is 

+ ( _ ( -rmu 2 + (k + h + k' L ) - Zt(uj) -k \ _1 , , 

° M " V -* -^2^ 2 + (fc + h + fcy - e+h J - (4 - i5j 

For this setup it again becomes difficult to evaluate explicitly the symplectic matrix S for 
the general case though it is straight-forward to evaluate it numerically to yield the steady 
state density matrix given by Eq. (2.28). 

We present some numerical results for this case. In our numerics we fix the following 
parameter values: m\ = 1.0, m 2 = 1.5, k — k± — k 2 — 1.0, k' L — k' R — 0.1. Moreover, we keep 
the temperature of the right reservoir fixed at T R = 1.0 and study the NESS for different values 
of AT = Tl — T r . We work in dimensionless units where K = k R = 1. The temperatures Tl, T r 
are of the order of the normal mode frequencies meaning indeed that the system operates in 
the quantum-mechanical regime. 
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The Green function in Eq. (4.15) is first obtained and then all the elements of the correlation 
matrix given by Eqs. (2.23) are evaluated by numerical integration. As examples we detail 
below the symmetric parts of the correlation matrices for the equilibrium case AT = and for 
AT = 4.0. 



Co = 



( 0.696 
0.294 




C, = 



\ 





1.851 
1.331 


-0.0294 



0.294 
0.670 





1.331 
2.241 
0.0196 







1.168 
- 0.0788 



0.0196 
2.491 
0.781 



/ 






- 0.0788 
1.67 

- 0.0294 N 



0.781 
4.558 / 



for AT = 



for AT = 4 



Note that the heat current across the chain is given by j = k{x\p2) /m,2 = 
{k/m-ijCu = —(k/mi)C23- For the above example we obtain j = 0.0196. 
As shown in Sec. (2.2) the NESS density matrix assumes the form: 



Ps = 



exp(—(p T Aip) 
4sinh(ai) sinh(a2) 



-k(x 2 pi)/m 1 = 



(4.16) 



where ip T = (xi,X2,pi,P2)- We next numerically determine ai,a 2 and the matrix A. To 
this end we need to construct the diagonal matrix D and the symplectic matrix S. The 
way of constructing these are described in Sec. (A): It requires the following four numerical 
procedures: 

(i) Find the eigenvalues and eigenfunctions of Cs- Then construct the matrix C l J 2 . 

(ii) Find the eigenvalues and eigenvectors of the matrix iC X J 2 J C^ 2 . There are four eigen- 
vectors which occur as complex conjugate pairs, (u^, wjf , wj, ui^), with corresponding eigen- 
values (—di, di, — d 2 , 0^2) ■ 

(iii) Evaluate the vectors vf = C^ 2 u^ , vf = Cg^uf and use Eqs. (A. 4, A. 20) to obtain 
the matrix V. The required symplectic transformation is then S = (JV) T . 

(iv) We evaluate a x = coth _1 (2rf 1 ), a 2 = coth~ 1 (2rf 2 ) and the steady state matrix 



Note that for AT = (T L 

Ps = Peg ~ e~ ms \ hence 



= S T Diag (01,02,01,02) <S>- 

Tr = 1) and for weak-coupling with reservoirs, we expect 



•-eq 



/ 


1 


0.5 








\ 




0.5 


1 


















0.5 







\ 











0.33. 


/ 



In Fig. (2) we depict the matrix elements A 33 , A 44 and A 34 as functions of the temperature 
difference AT. In the inset we have plotted the element Au and the heat current j = Cn/m 2 . 

The 2x2 diagonal blocks of the matrix A; i.e., An, A 12 , A 2 i, A 22 and A 33 , A 34 , A 43 , A 44 , 
can be obtained from the weak-coupling results in Sec. (3.2). First we obtain the normal 
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-M2 
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A 

A 33 


■ 


A 44 


♦ 


A 

A 34 




Figure 2: (color online). Plot of some relevant elements of the matrix A as a function of the 
temperature difference AT — T^ — Tr with constant Tr =1.0 while the other parameters are 
given in the text. The dashed lines depict results obtained from the analytical weak-coupling 
approximation. The inset shows both, the matrix element — A 14 and the linearly growing heat 
current j. 
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mode eigenvalues Q s and eigenfunctions Vi(s),s = 1,2, corresponding to the isolated system 
Hamiltonian The symplectic transformation is constructed by using Eq. (3.4). The matrix 
elements of D, given generally by Eqs. (3.9,3.10), takes the following form: 

i k'Wfis) , , hft s , i k' V 2 2 (s) , , m, . 

d s = = — L -^4 coth( — ) H 5 — coth( — ) , 

2^V 1 2 (s) + ^V 2 2 (s) y 2k B T L ^ 2k' L 2 V 1 2 (s) +k' R 2 V 2 2 (s) 2k B T R ' 

for s = 1,2. The weak-coupling results for A 33 , A 44 and A 34 have been plotted in Fig. (2) 
(dashed lines) and we detect an excellent agreement with the values obtained from precise 
numerics. 

5 Conclusions and outlook 

In summary, we have detailed the explicit construction of the reduced density matrix of 
the nonequilibrium steady states for two quantum transport problems, one involving non- 
interacting fermionic degrees of freedom and the other noninteracting bosonic degrees. The first 
setup concerns electron transport in a tight-binding lattice model composed of non-interacting 
electrons that are connected to non-interacting baths while our second setup focuses on heat 
transport across an arbitrary harmonic oscillator network connected to harmonic oscillator 
baths. For both these models the steady state correlations are known exactly from vari- 
ous approaches and are usually expressed in terms of nonequilibrium Green functions. We 
have demonstrated that for the Fermionic problem, the construction of the emerging time- 
independent steady state density matrix requires that one evaluates a particular unitary matrix 
while, likewise, for the Bosonic case, it requires finding an appropriate symplectic transforma- 
tion. 

For the limiting case of vanishingly weak coupling between intermediate system and reser- 
voirs, we show that the required unitary and symplectic transformations can be explicitly found 
and the resulting density matrices assume simple forms whose explicit expressions depend on 
the way the coupling strengths are made to vanish. For the case where the two baths possess 
the same temperatures (and chemical potentials for electron case) the weak coupling case yields 
a unique answer which is the expected equilibrium canonical (grand-canonical for electrons) 
distribution. This requires the assumption that the connecting reservoirs have sufficiently 
broad band- widths [16, 17]. 

The construction of the steady state density matrices required one to use "diagonal" rep- 
resentations [Eqs. (2.11,2.27)] and these are analogous to the eigenmode or normal mode 
representation of the Hamiltonian. In the equilibrium case and for weak coupling the den- 
sity matrix is ~ e~^ n and then the eigenmode representation is useful in the computation of 
equilibrium averages of various physical observables. Similarly, we expect that the "diagonal" 
representations of the nonequilibrium density matrix is as useful for computing averages in 
the NESS. Thus, for example, the Von Neumann entropy of the nonequilibrium steady state, 
defined as S = — Tr [ ps In p$ ] can be readily obtained from our findings. In particular one 
finds that: 

Sfcrmion = ~ Y.^li 1 ~ d s) hl(l - d a ) + d s kld s , . , 

Sboson = - Ef=i(4A - 1/2) Hdjh - 1/2) - (djh + 1/2) \n(d s /h + 1/2) , 1 • } 



20 



where {d s } are the "diagonalized" correlations defined via Eqs. (2.9, 2.25). 
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A Procedure to find the symplectic matrix S 

We here explain the general procedure to find the symplectic matrix S [24, 25]. We first 

ii 

consider the eigenvalue problem for the matrices iCgJCg and C$J ■ Note that the covariance 
matrix C$ is real- valued, symmetric and positive definite. Positive definitness is shown by 

y T Csy = y T Cy = {{ip T y) 2 ) ss > for arbitrary real column vector y. 

i i 

The matrix iCJJCJ is a Hermitian matrix. Therefore it possesses real eigenvalues as 
i i 

iCgJCgU = duo, where uj is the eigenvector. Taking the complex conjugate of both sides, we 

i i 

have the equation iCgJCgU* = —duo*. From this, if d is an eigenvalue, then — d is also an 
eigenvalue. 

Hence, we can start with the following equations 

ic\jc\u± = (A- 1 ) 

where u k are eigenvectors (u> k = oo k *) which have real eigenvalues ^fd k (d k > 0). These 
equations are equivalent to 



C s Jv± = ±id k v± (A.2) 



where the vectors v k are defined as 



v± = Cjut (A.3) 



We divide the vector vt into the real and imaginary parts as 



v 



v«±ivl (A.4) 



Then, Eq.(A.2) implies the two relations 

CsJvl = -v T k d k , (A.5) 

C s Jv{ = v*d k . (A.6) 

Because the matrix iCgJCg is Hermitian, we can normalize the vector oo k as 

(wf)t w ± = 2d~ k }8 k , k ,, (A.7) 

(<4M = °- ( A -8) 
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From (A. 3), the vector u k is expressed with vectors v k ' ! as 

wf = CsHv«±ivl). (A.9) 
Using this the Eqs.(A.7) and (A. 8) are written as 

{v?Tivl) T Cs*CP(v5±ivtt = {v£) T C?vg + (vl) T C?vl 

=F i MfC^v* - {v«) T C- s W k ] = 2d- k ,%, k , (A.10) 

(v^TivifC^C-Hv^Tivl,) = {v*) T C- s l v§-{vl) T C- s W k , 

=F i HfCsX' + K) T CsW kl ] = 0. (A.ll) 

From this, we find the following set of relations 

^) T C~ s l v§ = d-/<V, (A.12) 

(vlfCs'vl, = d-%, k ', (A.13) 

{v«) T C- s W k , = 0, (A.14) 

(vlfCsX' = 0. (A.15) 

Utilizing Eqs.(A.5) and (A. 6), the above relations can be recast as 

(vkfJ4 = hp, (A.16) 

(vi) T Jv§ = -6w, (A.17) 

(v*) T Jv§ = 0, (A.18) 

(vlfJvl, = 0. (A.19) 

We next define the 2N x 2N matrix V 

V = {v*,--- ,v§,v(,--- ,v T N ). (A.20) 

Using the matrix V, relations (A. 5) and (A. 6) can be simply written as 

C S JV = VJD, (A.21) 

where the matrix D is a 2N x 2N diagonal matrix 

D = Diag(di, • • • , d N , d 1 , ■ ■ ■ , d N ). (A.22) 

In addition, the relations (A.16)-(A.19) can be written with the matrix V as 

V T JV = J. (A.23) 

We now introduce the matrix S as 

S = (JV) T (A.24) 
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One can prove that the matrix S satisfies the symplectic relation, namely: 



SJS T = V T J T JJV 

= -V T J T V = V T JV = J, (A.25) 

SC S S T = V T J T C S JV 

= V T J T VJD 

= -V T JVJD 

= -J 2 D = D, (A.26) 

where we used Eqs.(A.21) and (A. 23). 

To evaluate the symplectic matrix S numerically, we first solve eigenvalue problem (A.l) 
to obtain the eigenfunction u^. Next, we normalize them as in (A. 7), and find v^' 1 . Finally, 
constructing the matrix V as in (A. 20), One obtains the symplectic matrix as in (A. 24). 
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